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ABSTRACT 


Progressive failure analysis (PFA) of non-crimp fabric (NCF) composite laminates 
subjected to low velocity impact loads was performed using the Complete STress 
Reduction (COSTR) damage model implemented through VUMAT and UMAT41 
user subroutines in the frame works of the commercial finite element programs 
ABAQUS/Explicit and LS-DYNA, respectively. To validate the model, low velocity 
experiments were conducted and detailed correlations between the predictions and 
measurements for both intra-laminar and inter-laminar failures were made. The 
developed material and damage model predicts the peak impact load and duration very 
close with the experimental results. Also, the simulation results of delamination 
damage between the ply interfaces, in-plane matrix damages and fiber damages were 
all in good agreement with the measurements from the non-destructive evaluation 
data. 
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INTRODUCTION 


Carbon fiber composite materials are being widely used in aerospace applications 
and are being slowly introduced into the automobile industry to achieve better fuel 
economy, increased durability and improved driving experience. Optimizing the 
weight of an automobile is a significant factor in achieving higher fuel economy, and 
to that end, composite materials are very promising due to their low densities. Among 
the carbon fiber composite materials under consideration, the non-crimp fabric (NCF) 
preform is an attractive alternative to the traditional pre-pregnated tapes due to its low 
manufacturing cost and high efficiency. The NCF material consists of multiple layers 
of unidirectional plies oriented at varying angles which are held together by transverse 
stitches. A typical carbon fiber NCF composite laminate is shown in Figure 1. 



Figure 1 : Photograph of a typical Non-Crimp Fabric (NCF) laminate. 

Since the use of NCF composites is gaining momentum, researchers are 
developing novel computational methods to evaluate the performance of these 
materials under varying load conditions. Researchers in Refs. [1-3] have developed 
micro, macro and meso-scale material models to evaluate some of the material 
properties and strengths, which are extremely difficult to obtain by standard material 
characterization tests. In Refs. [4-5], researchers have developed a meso-scale plane 
stress damage model to simulate damage in NCF cross-ply laminates subjected to 
tension and compressive loads. The influence of out-of-plane imperfections in the 
NCF laminate due to the misalignment of the fibers was investigated in Ref. [5], The 
extent of damage in a NCF composite sandwich panel under impact load was 
established using an ultrasonic C-scan in Ref. [6], 

Even though a good amount of research has been performed in simulating the 
damage progression in tape laminates subjected to impact loads [7-11], less effort has 
been dedicated to NCF composite laminates. Hence, in the present investigation, a 2D 
stress damage model called Complete STress Reduction (COSTR), previously 
developed for fiber and matrix damage predictions in tape laminates by the authors 
from NASA [12] for in-plane loads, was extended to out-of-plane loads for thin NCF 
laminates. This damage model will be implemented in the framework of the 


commercial software packages ABAQUS and LS-DYNA. Though this approach was 
first developed for use in Abaqus for the aerospace industry, because the automotive 
industry relies heavily on LS-DYNA, the damage models developed for NCF 
composites will be incorporated into both codes. In this manner it is anticipated that 
large real world problems such as the progressive failure of aircraft wings and 
fuselages, and crash of automobiles, can be simulated. NCF material preparation and 
low velocity dynatup impact tests were conducted at the General Motors R & D 
Facility. Progressive failure analysis (PFA) of the NCF laminate, impacted by a 
hemispherical impactor at 5J, 10J & 12J energy levels was performed and analysis 
predictions for the overall panel response, such as impact load versus time duration, 
were compared with the experimental results. Fiber damages within the plies and 
delamination damages between the layers were compared with X-ray and CT-scan 
data. 


PROGRESSIVE FAILURE ANALYSIS - DAMAGE MODELS 

During the low velocity impacts, composites are subjected to both intra-laminar 
and inter-laminar failures. Intra-laminar failures were studied previously for in-plane 
loads utilizing the COSTR damage model developed and implemented through a 
VUMAT user subroutine in ABAQUS/Explicit [12] for tape laminates. The inter- 
laminar failure, referred to here as delamination failure occurring between the layers 
of a composite, poses a major challenge for damage prediction and thus requires an 
advanced modeling strategy. A brief description of the damage model will be outlined 
here. For a detailed description of the model, readers are referred to refs. [12-13]. 

Intra-Laminar Damage Model 

The COSTR damage model for intra-laminar damage utilizes the Hashin-Rotem 
unidirectional failure criteria [14] to detect fiber and matrix damages in a lamina. The 
failure criteria are expressed in terms of the in-plane stresses o ip the strengths in 
tension, X, compression, Y, and shear, Sxy, respectively. Two damage indices d F and 
d M are used to represent the fiber and matrix damage modes, respectively. Each 
damage index has the value of zero (no damage) or 1.0 (complete damage). The 
previous stress degradation model, adopted to degrade the axial (an), transverse (022), 
and shear ( 012 ) stresses as a function of axial (£ 11 ), transverse (£22), and shear (£ 12 ) 
strains respectively, after they reach their failure limits at a material point is shown in 
Figure 2(a). When the damage parameter d M attains a value of 1.0, both the transverse 
and the shear stresses were reduced to zero instantaneously. Similarly, when the 
damage parameter d F attains a value of 1.0, all of the in-plane and shear stresses at that 
material point were made zero. This type of instantaneous stress degradation to zero is 
characterized as a single-step or instantaneous degradation approach and hence the 
model was assigned the name Complete STress Reduction (COSTR) damage model. 
This type of stress degradation approach is suitable for performing PFA of un- 
configured laminated composite panels under quasi static loading conditions. But, 
while analyzing configured panels (skin-stringer) or structures under dynamic loading 
conditions, this stress degradation approach generates an extremely large energy 
release in the model thereby causing all the elements along the crack paths to fail 



instantaneously. This phenomenon has been studied and reported in Ref. [14]. To 
circumvent the above difficulties, the fiber stresses in the failed material point are 
reduced to a nominal percentage of the original strength (20% or so), and these 
stresses are maintained until the strain in the fiber direction reaches a value of 2.5 
times the elastic strain (£ el n) of the fiber. However, the transverse and shear stresses 
are still degraded to zero instantaneously after the failure as they do not pose 
difficulties in the simulation. The stress degradation model used in the present study is 
shown in Figure 2(b). 




a) Previous Approach b) Current Approach 

Figure 2. Stresses degradation approaches in COSTR damage model. 

Another characteristic of the COSTR damage model that differs from the standard 
ABAQUS/LS-DYNA damage model is the elimination of the influence of Poisson’s 
ratio in evaluating the axial stress in the fibers when the matrix has failed due to the 
transverse tension/compression or shear stress. It has been observed during the study 
that in some plies of a laminate, the transverse strain, shear strain or the combination 
of transverse and shear strains could result in significant distortion of the element 
without significant stress in the fiber direction. An example of this phenomenon is a 
laminate consisting of 90° plies and loaded in the direction perpendicular to the fibers. 
To avoid excessive distortion of an element, a quadratic strain criterion, was 
formulated and implemented in the damage evolution procedure of the PFA 
methodology to eliminate badly distorted elements from the simulation. The matrix 
strain failure criterion (MSFC) is defined as 


MSFC = a /(^ 22) 2 + (r l2 ) 2 

When the MSFC at a material point reaches a significant value, approximately 0.25, 
the material point is deleted. Once a material point is deleted it cannot be reactivated. 
The ABAQUS and LS-DYNA explicit solvers check for elements where all of the 
material points have been flagged as deleted and then, these elements are removed 
from further computations, thereby simulating virtual cracks in the finite element 
model. The COSTR damage model has been implemented through the VUMAT and 
UMAT41 [15] user written subroutines in the ABAQUS and LS-DYNA codes, 
respectively. 


Inter-Laminar Damage Model 


The simulation of inter-laminar damage, i.e., delamination, between the 
sublaminate layers was accomplished using the cohesive zone model available in the 
ABAQUS and LS-DYNA computer programs. The onset of delamination was 
determined based on the inter-laminar quadratic normal and shear stress criterion and 
the delamination growth was based on a critical fracture energy criterion. Damage 
was modeled as an irreversible process by including a damage parameter and this 
parameter is directly related to the dissipated fracture energy. A detailed explanation 
of the cohesive zone model can be found in Ref. [13]. Among the analysis parameters 
for the cohesive layers, stiffness of the cohesive layer is an important parameter for the 
simulation, as lower stiffness values can make the laminate too compliant and higher 
stiffness values can cause spurious oscillations of the tractions in the element 
introducing numerical difficulties. Ref [16] provides the guidelines for stiffness 
determination and one such guideline for the mode-1 direction is given below: 


where E 3 is the Young’s modulus of the laminate in the thickness direction, t is the 
larger of the sub-laminate thicknesses above or below the cohesive layer, and a is a 
parameter that is much larger than 1. A value of a equal to 50 is recommended in Ref. 
[15] and was used in the current study to obtain a stiffness of the cohesive layer. In 
calculating the stiffnesses K u and K m for the mode-1 and mode-2 directions, E 3 is 
replaced in the above equations with the shear moduli Gn and G13 of the laminate, 
respectively. 


DESCRIPTION OF NCF COUPON AND TEST SETUP 

An eight layer carbon fiber non-crimp fabric with a quasi-isotropic layup (0/ 
-45/45/90/90/45/-45/0) from Sigmatex was used to manufacture plates with 
dimensions 431.8x431.8x 1.8 mm at the General Motors Research Labs using resin 
transfer molding. The NCF fabric was made from T700 grade fiber. An epoxy resin 
with a fast curing agent was used for the resin transfer molding of the composite. The 
volume fraction of carbon fiber was determined to be close to 48%. Basic material 
characterization tests for stiffness and strength were conducted following ASTM 
standards to obtain the material properties in tension, compression, and shear. The 
material properties are listed in Table I. These material properties were validated by 
correlating the results with a three point bend test. For the low velocity impact tesing, 
square specimens of 101.6 mm wide were cut from the manufactured plates of NCF 
and were subjected to an impact using an INSTRON dynatup testing machine, 
equipped with a digital acquisition system. The test-specimen was restrained using a 
clamping fixture (upper and lower plates) which exposed a circular area of 76.2 mm 
diameter, and then impacted by a rigid tup at the center of the exposed area. A 
hemispherical steel impactor with a diameter of 25.4 mm and a weight of 7.0 kg was 
used. Three levels of energy input (5J, 10 J and 12J) were studied. The test setup is 
shown in Figure 3. 



TABLE I. MATERIAL PROPERTIES OF NCF COUPONS 


Properties 

T700/Epoxy 

Description 

En (GPa) 

102.0 

Young’s modulus in fiber direction 

E 22 (GPa) 

7.08 

Young’s modulus in the transverse direction 

G i2 (GPa) 

3.16 

In-plane shear modulus 

S„ (GPa) 

1.63 

Tensile strength in the fiber direction 

C„ (GPa) 

1.43 

Compressive strength in the fiber direction 

S 22 (GPa) 

0.068 

Tensile strength in the transverse direction 

C 22 (GPa) 

0.30 

Compressive strength in the transverse direction 

S 12 (GPa) 

0.05 

In-plane shear strength 

v i 2 

0.4 

Poisson’s Ratio 

G lc (J/m2) 

1062.0 

Critical Fracture toughness in Mode 1 Direction 

Guc Gnic (J/m2) 

1276.0 

Critical Fracture toughness in Mode 2 & 3 directions 


The cohesive damage model available in ABAQUS and LS-DYNA was used in 
simulating delaminations between the sublaminate layers. The interfacial stiffness and 
strength parameters used in the analysis are given in Table II. 


TABLE II. INTERLAMINAR PROPERTIES AND STRENGTH DATA 


Properties 

T700/Epoxy 

Description 

K„ (GPa) 

781.4 

Interfacial Stiffness in Mode 1 direction 

Ki 2 (GPa) 

360.1 

Interfacial Stiffness in Mode 2 direction 

K 22 (GPa) 

360.1 

Interfacial Stiffness in Mode 3 direction 

S33 (GPa) 

0.06895 

Interfacial strength in Mode 1 direction 

S 13 (GPa) 

0.08274 

Interfacial strength in Mode 2 direction 

S 32 (GPa) 

0.08274 

Interfacial strength in Mode 2 direction 


Dyna Tup 



Figure 3: Assembly of Non-Crimp Fabric (NCF) laminate & AL plate with cutouts. 


FINITE ELEMENT MODEL DESCRIPTION 


A three-dimensional (3D) finite element model of the square plate was developed 
to simulate the coupled inter- and intra-laminar damage. The thickness of the model 
was divided into five sub-laminate sections with a cohesive layer of zero thickness 
between each of them. The pictorial representation of the plate model, which consists 
of sub-laminate sections and the cohesive layers, is shown in Figure 4(a). The finite 
element mesh for each of the sub-laminates and the cohesive layers is shown in Figure 
4(b). Each sublaminate layer was modeled using one solid element in the thickness 
direction. The sub-laminate sections were discretized using the 8-node continuum 
shell reduced integration element (SC8R) available in the ABAQUS™ element library 
[13] and continuum solid elements in LS-DYNA. The SC8R element is a 3D element 
with translational degrees-of-freedom (DOF) using only linear interpolation functions. 
It differs from a 3D solid element in that it employs composite lamination theory and 
therefore provides an ability to discretize 3D bodies with multiple sublaminates 
through the thickness. The element size used in discretizing the planar regions of the 
sub-laminates ranges from 0.5 mm by 0.5 mm in the center under the impactor to a 
size of 1.5 mm by 1.5mm close to the edges of the plate. The cohesive layers were 
modeled using a zero thickness ABAQUS™ 8 -node cohesive element known as 
COH3D8 [16], and in LS-DYNA, cohesive elements were defined with material 
model MAT_COHESTVE_MIXED_MODE. An in-plane discretization of 0.5 by 0.5 
mm was used in the cohesive layers. The sub-laminate sections and the cohesive 
layers were connected using ABAQUS™ and LS-DYNA “tie” boundary conditions, 
in each approach respectively, to maintain the continuity between the sub-laminate 
sections before the damage occurs. The assemblage of the five sub-laminate sections 
and four cohesive layers through the thickness defines the complete plate model. 



Sub-Laminate Layer 





Cohesive 

Layer 


(a) Schematic Diagram of Sub-laminate & 
Cohesive layers 


(b) Finite Element Mesh 


Figure 4: Modeling procedure of NCF plate & finite element mesh. 

The nodes on the sublaminate surfaces which are in contact with the top and 
bottom support plates were constrained in the out of plane direction only. Additional 
boundary conditions were applied to the nodes along the vertical and horizontal center 


lines of all the sublaminate layers as shown in Figure 5 in order to suppress spurious 
element distortions due to shock waves. The plate was impacted with an impactor at 
an initial velocities of 1.12 m/sec, 1.68 m/sec and 1.85 m/sec to produce the desired 
5J, 10J and 12J energy levels. 



Figure 5: Finite element mesh with prescribed boundary conditions. 

The impact simulation was performed using the ABAQUS™ and LS-DYNA 
explicit double precision solvers for each of the energy levels. Since the plate is 
severely constrained and the effective width of plate was rather small, the elements 
under the impactor can show spurious hour-glass modes due to the under integrated 
elements. Hence, an additional viscous stiffness was added in all the sublaminate 
elements to eliminate zero energy modes in the plate in the respective codes. 

RESULTS AND DISCUSSION 

PFA results such as impact load history and delamination damage at certain 
critical interfaces of the NCF plate for the impact energy levels of 5J, 10J & 12J are 
compared with experimental test results in the following sections. For the sake of 
brevity, the comparative simulation results for LS-DYNA and the ABAQUS/COSTR 
damage model were presented only for the 5J test case. 

5J Impact Energy 

The load history curves obtained from the PFA and two tests for the 5J energy 
level impact cases are presented in Figure 6. From Figure 6, it can be noticed that the 
contact duration and the peak impact load obtained from the analysis (green) are in 
good agreement with the test data (purple). The raw load history curve from the 
analysis was oscillatory, due to the contact disturbances occurring in the simulation. 
These oscillations in the load were smoothened by passing it through filters provided 
in the ABAQUS post processing tool. 




Figure 6: Comparison of load versus time curve from ABAQUS/VUMAT simulation and dynatup 

experiment 

Figure 7 presents the comparison of load history curves for the LS-DYNA user 
defined COSTR damage model and test results. It can be observed that the damage 
model implemented in LS-DYNA provides excellent correlation with the 
experimental results and also agrees well with the ABAQUS/COSTR simulation. A 7 
point averaging available in the LSPOST program was used to smooth the contact 
force curves. 



Figure 7 : Comparison of load versus time curve from LS-DYNA UMAT 41 simulation and 

dynatup experiment 


The simulated delamination at certain critical ply interfaces for the 5J impact load 
is presented in Figure 8, along with the NDE ultrasound measurements. The general 
orientation of the delamination damage showed fair agreement at the various 
interfaces in comparison with the NDE data. The layers in the NCF plate are 
numbered such that the layer that comes in contact with the impactor was numbered as 
layer 1 and the bottom most layer was numbered as layer 8. The cummulative 
delamination through the entire thickness from the ABAQUS/VUMAT simulation 
was compared with the damage measured using ultrasound NDE techniques and is 
shown in Figure 9(a). From the figure, good agreement is observed between the 
predictions and the experimental measurements. Similarly, the simulation results for 
the cumulative delamination in the laminate from the LS-DYNA/UMAT41 simulation 
was compared with ultrasound NDE results and presented in Figure 9(b). A good 
correlation with the experiment and both the simulation results was observed, 
confirming a successful implementation of the COSTR damage model in the LS- 
DYNA code. In all of the delamination pictures below, the region in white indicates a 
complete delamination of interfaces, the region in red indicates areas which are 95% 
of complete delamination and the region in blue implies an intact bonding of the 
interfaces. 

No fiber failure was predicted in any of the plies in the simulations, matching the 
experimental observations. The matrix damage predicted from the simulation in the 
elements was not significant, and hence, is not presented. 
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Figure 8: Delamination at certain interfaces. Left hand side - PFA, Right hand side - NDE technique. 









(a): Left hand side - AB AQU S/VUM AT simulation, Right hand side - NDE technique. 
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(b): Left hand side - LS-DYNA /UMAT41 simulation, Right hand side - NDE measurements. 
Figure 9: Overall delamination of NCF laminate for a 5J impact load. 


10J Impact Energy 

The load history curves and delamination damage obtained from the 
ABAQUS/COSTR PFA were compared with the four experiments with 10J impact 
energy and are presented in Figure 10. 





Figure 10: Load history curve for 10J impact energy case. 

It can be seen from the figure that the contact duration and the impact load 
obtained from the analysis are in very good agreement with the test data. Similar to the 
5J impact energy case, the load history curve obtained from the analysis was 
oscillatory, and by passing through the multi-pass filters provided in the ABAQUS 
post processing tool, these oscillations were smoothed out. 

The delamination at certain critical ply interfaces for the 10J impact load is 
presented in Figure 11. The ultrasound scans for this specimen were taken looking up 
from the bottom, starting from the 8 th layer (non-impact side). The general orientation 
of the delamination at the interface of layers 6-7 was predicted fairly well. However 
the NDE data at other interfaces show poor correlation with the simulation, but this 
may be due to the masking effect of the damage at the interface of layers 6-7. The 
areas of cumulative delamination of the laminate, determined from PFA, and the 
experimental measurements from the NDE ultrasound are presented in Figure 12. 
These show a good agreement and are a better assessment of the total damage area in 
the laminate. 
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Figure 11: Delamination at certain interfaces. Left hand side - PFA, Right hand side - NDE 

technique. 






Figure 12: Overall delamination of NCF laminate for a 10J impact load. Left hand side - PFA, 

Right hand side - NDE technique. 

For the 10J impact energy level, a significant amount of matrix damage was 
predicted in all of the plies as presented in Figure 13. The region with matrix damage 
is indicated in red and the undamaged regions are indicated in blue. The matrix 
damage predicted in the plies is due to the in-plane transverse tension/compression 
stress, or shear stresses, reaching the failure limit either individually or in combination. 
From Figure 13, it is seen that the matrix damage was occurring parallel to the fibers 
in the 0° and 90° plies (see the layer-8 and layer-4&5 plots). However for the 45° plies, 
the damage was oriented parallel to, or and at some angle to the fibers (see the layer-6 
and layer-3 plots). No fiber damage was predicted in any of the plies in the simulation 
and none was observed in the experiments. 
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Figure 13: Matrix damage in plies - PFA. 








12J Impact Energy 


The load history curves obtained from PFA and four experiments for the impact 
test case of 12J energy level are presented in Figure 14. The magnitude of the impact 
load and the contact duration from the simulations are in good agreement with 



Figure 14: Load history curve for 12J impact energy case. 


the test data. The NCF plate had extensive matrix damage in all of the plies and fiber 
failure in some plies. The test data shown in Figure 13 indicates that this damage 
occurs at around 3.1 milliseconds, accompanied by a load drop. The PFA simulation 
predicts similar behavior in the load history curve. 

The delamination at certain critical ply interfaces for the 12J impact load is 
presented in Figure 15. As the non-destructive evaluation results for the delamination 
damage at each of the interfaces was not available, only simulation results are 
presented. The cumulative delamination damage in the NCF plate predicted from the 
PFA was compared with the NDE results and is presented in Figure 16. From this 
figure, we can see a good agreement in the overall size of the delamination between 
the PFA and experimental results. The length of delamination damage from the test 
data (69 mm) was larger than the predicted delamination from the PFA data. This is 
due to the fact that few fiber tows in the outer layer of the laminate from the non- 
impact side debonded from the rest of the laminate and the PFA was unable to capture 
this phenomenon as layers 7 & 8 were modeled with one element in the thickness 
direction thereby restricting the relative displacement between them. However, the 
delamination damage in the rest of the plies matches well with the PFA data. Again, 
as mentioned in the previous sections, the white region indicates a complete 
delamination of the interfaces, the red region indicates areas which are 95% of 
complete delamination of the interface, and the blue region denotes intact bonding of 
the interfaces. 
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Figure 15: Delamination at certain interfaces. Left hand side - PFA, Right hand side - NDE 

technique. 



Figure 16: Overall delamination of NCF laminate for a 12J impact load. Left hand side - PFA, 

Right hand side - NDE technique. 

For the 12J impact energy level, a significant amount of matrix damage was 
predicted from the PFA in all the plies, and is presented in Figure 17. The region with 
matrix damage is indicated in red while the undamaged regions are indicated in blue. 
It can be observed from Figure 17 that the matrix damage is again occurring parallel to 
the fibers in the 0° and 90° plies (see layer-8, layer-4&5). However, in the 45° plies 





(layer-2 and layer-3) , the damage is oriented in parallel to, or at some angle relative to 
the fibers. Note the circumferential damage in layer- 1 at the circular clamping 
boundary. 



O 

! h. Figure 17: Matrix damage in plies - PFA. 


Figure 18 presents the comparison of fiber damage in the outer 5 plies (from layer- 
4 to layer-8 where layer- 8 is the bottom most layer) predicted from simulation and CT 
scan measurements. Complete fiber failure regions are indicated in white and partial 
fiber failure areas are indicted in blue. Undamaged areas are indicated in red. Note 
from the figure that the fiber damage predicted in the simulation (white regions) 
matches well with the fiber failure determined from the CT scan data. The fiber 
damage in layers 7 & 8 are identical in shape and size in the simulation due to the fact 
that the relative displacements between these two layers have been restricted on 
account of modeling procedure adopted in this work. 
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Figure 18: Fiber damage in plies - PFA. 


CONCLUSION 


Low velocity impact tests and progressive failure analysis of thin NCF 
laminated composites were successfully performed for 5J, 10J and 12J energy levels. 
The Complete Stress Reduction (COSTR) damage model was newly implemented in 
the LS-DYNA framework using a user subroutine based on the previous 
ABAQUS/VUMAT implementation. Both of the implementations were assessed in 
predicting the fiber, matrix, and delamination damages under the low velocity 



transverse impact and found to be in good agreement with the damage identified 
through NDE technique used during testing. 
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